##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)
library(BSgenome.Hsapiens.UCSC.hg38)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/all/testis_combined_peak.combineRNA.qc.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac_v3/all"

}


###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

addArchRThreads(threads = 10)
addArchRGenome("hg38") 

##########################################################################################
 
a <- load(comine_data_file)

##########################################################################################
# 调用macs2 
pathToMacs2 <- findMacs2()
#pathToMacs2 <- "~/miniconda3/envs/archr/bin/macs2"

# call peak
testis_combined_peak_combineRNA <- addGroupCoverages(testis_combined_peak_combineRNA,groupBy = "cell_type",force = TRUE)  ## 创建伪批量重复
testis_combined_peak_combineRNA <- addReproduciblePeakSet(testis_combined_peak_combineRNA,groupBy = "cell_type",force = TRUE)  ## 迭代重叠峰
testis_combined_peak_combineRNA <- addPeakMatrix(testis_combined_peak_combineRNA,force = TRUE)

## motif
testis_combined_peak_combineRNA <- addMotifAnnotations(ArchRProj = testis_combined_peak_combineRNA, motifSet = "cisbp", name = "Motif",force = TRUE)
testis_combined_peak_combineRNA <- addBgdPeaks(testis_combined_peak_combineRNA,force = TRUE)
testis_combined_peak_combineRNA <- addDeviationsMatrix(
  ArchRProj = testis_combined_peak_combineRNA, 
  peakAnnotation = "Motif",
  force = TRUE
)

## peak2gene
testis_combined_peak_combineRNA <- addPeak2GeneLinks(ArchRProj = testis_combined_peak_combineRNA , useMatrix = "GeneExpressionMatrix")

## testis_combined_peak_combineRNA
testis_combined_peak_combineRNA@cellColData$seurat_clusters <- as.character(testis_combined_peak_combineRNA@cellColData$seurat_clusters)
testis_combined_peak_combineRNA@cellColData$seurat_clusters <- ifelse( 
    testis_combined_peak_combineRNA@cellColData$seurat_clusters %in% c(12,13) , 17 ,
    testis_combined_peak_combineRNA@cellColData$seurat_clusters
    )
testis_combined_peak_combineRNA@cellColData$seurat_clusters <- as.numeric(testis_combined_peak_combineRNA@cellColData$seurat_clusters )
testis_combined_peak_combineRNA@cellColData$seurat_clusters <- factor(testis_combined_peak_combineRNA@cellColData$seurat_clusters)


## 输出
out_file <- paste0( out_path , "/testis_combined_peak.combineRNA.qc.Rdata" ) 
save( testis_combined_peak_combineRNA, file = out_file )
saveArchRProject( testis_combined_peak_combineRNA , out_path )

##########################################################################################
## peak的counts矩阵输出为seruat对象
all_peak <- data.frame(testis_combined_peak_combineRNA@peakSet)
all_peak$cell_type <- names(testis_combined_peak_combineRNA@peakSet)
cell_id <- testis_combined_peak_combineRNA@cellColData[,c( "cell_type" , "seurat_clusters")]
all_peak <- merge( all_peak , unique(cell_id[ , c("cell_type" , "seurat_clusters")]) , by = "cell_type")
all_peak <- all_peak[, c(2:ncol(all_peak) , 1) ]

out_file <- paste0( out_path , "/testis_combined_peak.qc.tsv" ) 
write.table( all_peak , out_file , row.names = F , quote = F , sep = "\t" )
